/* gcc -O2 parse.c -o parse.x -lm */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define max(A,B) ((A) > (B)) ? (A) : (B)
#define MAXBINS 25


/* void moments(double *x,int n,double nv){ */
/*   x[1] /= x[0]; */
/*   x[2] /= x[0]; */
/*   x[2] = sqrt((x[2]-x[1]*x[1])/x[0]); */
/*   x[1] /= nv; */
/*   x[2] /= nv; */
/*   if(n == 0){ */
/*     printf(" %10.7e %5.1e ",x[1],x[2]); */
/*   }else{ */
/*     printf(" %5.3e %5.1e ",x[1],x[2]); */
/*   } */
/*   x[0] = x[1] = x[2] = 0; */
/* } */


void initparam(int *nvmin,int *nvmax,int *minseed,int *maxseed,float *taumin,float *taumax,float *dtau){
    FILE *fp=fopen("parameters","r");
    fscanf(fp,"%d\t%d\t%d\t%d\t%f\t%f\t%f",nvmin,nvmax,minseed,maxseed,taumin,taumax,dtau);
    fclose(fp);
}


/*
 inputs an integer n and converts it into a sequence of characters in *str,
 starting at str[i] and returning the integer i' where "n" ends in *str
 (str[i]='\0').
 */
int inttostring(int n,int i,char *str){
    if(n < 0){
        str[i++] = '-';
        i = inttostring(-n,i,str);
    }else{
        if((n/10) > 0) i = inttostring(n/10,i,str);
        str[i++] = '0' + n%10;
        str[i] = '\0';
    }
    return(i);
}

FILE *opennewfile(int nv,float tau,int seed){
    int k,i=0;
    FILE *fp;
    char fname[100];
    char s;
    int runpow;
    float runfac;
    
    i = inttostring(nv,i,fname);
    fname[i++] = '/';
    fname[i++] = 'V';
    fname[i++] = 'E';
    fname[i++] = 'O';
    fname[i++] = 'g';
    fname[i++] = 'p';
    fname[i++] = 'u';
    fname[i++] = '1';
    fname[i++] = '3';
    fname[i++] = '_';
    i = inttostring((int)(100*tau+0.01),i,fname);
    fname[i++] = '_';
    i = inttostring(seed,i,fname);
    fp = fopen(fname,"r");
    if(fp != NULL){
        k = 2;
        do{
            fscanf(fp,"%c",&s);
            k -= (s == '\n');
        } while( k );
        //            printf("#  %s %c\n",fname,s);
    }
    return(fp);
}

FILE *openoutfile(int nv){
    int k,i=0;
    FILE *fp;
    char fname[100];
    
    i = inttostring(nv,i,fname);
    fname[i++] = '/';
    fname[i++] = 'P';
    fname[i++] = 'D';
    fname[i++] = 'F';
    fname[i++] = '3';
    fname[i++] = 'S';
    fname[i++] = 'p';
    fname[i++] = 'i';
    fname[i++] = 'n';
    i = inttostring(nv,i,fname);
    fp = fopen(fname,"w");
    return(fp);
}




int main(void){
    FILE *fp,*fp_out;
    int nv,nnv,nvmin,nvmax;
    int i,inst,inst_old,minseed,maxseed,seed1,run,magnit,magnit_old;
    float x,taumin,taumax,dtau,tau,mintime;
    int distr[MAXBINS],nvflag;
    double minenergy;
    double bestenergy,lowE,highE,dE,sigma,avertime[3],ener[3],allenergy;
    
            for(i=0; i<MAXBINS; i++) distr[i]=0;
    initparam(&nvmin,&nvmax,&minseed,&maxseed,&taumin,&taumax,&dtau);
    for(nv=nvmax; nv>=nvmin; nv--){
        nvflag = 0;
        for(i=0; i<3; i++) ener[i] = 0;
        for(i=0; i<3; i++) avertime[i] = 0;
        for(tau=taumin; tau<=taumax; tau+=dtau){
            for(seed1=minseed; seed1<=maxseed; seed1++){
                fp = opennewfile(nv,tau,seed1);
                if(fp != NULL){
                    nvflag = 1;
                    //      40       9323 -690.0000000000            2898        33167     2   14.0000000000       0
                    while(fscanf(fp,"%d%d%lf%g%*d%*d%*lf%d",&nnv,&inst,&minenergy,&mintime,&magnit) != EOF){
                        //                    printf("%d %d %d %lf %g %d\n",seed1,nnv,inst,minenergy,mintime,magnit);
                        bestenergy = (double)minenergy/nnv/nnv;
                        ener[0] += 1;
                        ener[1] += bestenergy;
                        ener[2] += pow(bestenergy,2);
                        avertime[0] += 1;
                        avertime[1] += mintime;
                        avertime[2] += mintime*mintime;
                    }
                    fclose(fp);
                }//if fp
            }//for seed
        }// for tau
        //              printf("count=%f %f %f\n",ener[0],ener[1],ener[2]);
        if(nvflag){
            fp_out = openoutfile(nv);
            avertime[1] /= ener[0];
            ener[1] /= ener[0];
            ener[2] /= ener[0];
            sigma = sqrt(ener[2]-ener[1]*ener[1]);
            dE = sigma*16/MAXBINS;//measure only 8 sigma wide
            fprintf(fp_out,"# %5d %8.0f %15.10e %15.10e %15.5e %15.5e \n",nv,ener[0],ener[1],sigma/sqrt(ener[0]),avertime[1],dE);
            printf("%5d %8.0f %15.10e %15.10e %15.5e %15.5e \n",nv,ener[0],ener[1],sigma/sqrt(ener[0]),avertime[1],dE);
//            for(i=0; i<MAXBINS; i++) distr[i]=0;
            
            for(tau=taumin; tau<=taumax; tau+=dtau){
                for(seed1=minseed; seed1<=maxseed; seed1++){
                    fp = opennewfile(nv,tau,seed1);
                    if(fp != NULL){
                        while(fscanf(fp,"%d%d%lf%g%*d%*d%*lf%d",&nnv,&inst,&minenergy,&mintime,&magnit) != EOF){
                            bestenergy = (double)minenergy/nnv/nnv;
                            //          minenergy = (allenergy+minenergy)/2;
                            i = (int)((bestenergy-ener[1])/dE+MAXBINS/2);
                            if((i > 0)&&(i<MAXBINS)){
                                ++distr[i];
                                ++distr[0];
                            }
                        }
                        fclose(fp);
                    }//if fp
                }//for seed
            }// for tau
            
            for(i=1; i<MAXBINS; i++){
                if(distr[i] != 0) {
                    x = (i-MAXBINS/2+0.5)*dE/sigma;
                    fprintf(fp_out,"%20.10e   %20.10e   %d \n",x,distr[i]*sigma/dE/distr[0],distr[i]);
                }
            }//if MAXBINS
            fclose(fp_out);
        }// if nvflag
    }//for nv
    return(0);
}

/********** END of MAIN           *** *****************************/
